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ABSTRACT 

We present a comprehensive numerical study of the dynamics of relativistic axisymmetric 
accretion tori with a power-law distribution of specific angular momentum orbiting in the 
background spacetime of a Kerr black hole. By combining general relativistic hydrodynamics 
simulations with a linear perturbative approach we investigate the main dynamical properties 
of these objects over a large parameter space. The astrophysical implications of our results 
extend and improve two interesting results that have been recently reported in the literature. 
Firstly, the induced quasi-periodic variation of the mass quadrupole moment makes relativis- 
tic tori of nuclear matter densities, as those formed during the last stages of binary neutron 
star mergers, promising sources of gravitational radiation, potentially detectable by interfero- 
metric instruments. Secondly, p-mode oscillations in relativistic tori of low rest-mass densities 
could be used to explain high frequency quasi-periodic oscillations observed in X-ray binaries 
containing a black hole candidate under conditions more generic than those considered so far. 

Key words: accretion discs - general relativity - hydrodynamics - oscillations - gravitational 
waves 



1 INTRODUCTION 

Relativistic tori (i.e. geometrically thick discs) orbiting around 
black holes are expected to form in a number of different scenar- 
ios, such as after the gravitational collapse of the core of a ro- 
tating massive star (M > 25Mq) or after a neutron star binary 
merger. State-of-the-art nu merical simulations of these scen arios, 
both in Newtonian physics fcxell200ll iRuffert & Jankdl200ll) . as 
well as in a relativistic framework ( Shibata et al. 2003), have shown 
that under certain conditions, a massive disc may be produced. 
This disc surrounds the newly formed black hole, has large aver- 
age rest-mass densities and highly super-Eddington mass fluxes, 
and its angular momentum obeys a sub-Keplerian distribution. Ac- 
cretion tori of much smaller accretion rates have also been con- 
sidered in general relativistic magnetohydrodynamic simulations 
of accretion flows onto Kerr black holes jbe Villiers et al] |2003|: 
i Gammie. McKinnev & T6thl2003l:lGammie. Shapiro & McKinnev 

There are several reasons for considering these objects astro- 
physically interesting. One of the most important reason has to do 
with the idea that an accreting torus around a black hole gener- 
ated after neutron star coalescence could be the progenitor source 
for short duration, (i.e. < 2 s) gamma-ray bursts (GRBs), where 
the relativistic fireball expands along the rotati on axis with very 
low baryonic contamination iGoodmarl 1 986llPaczvnskil 1 98d) (see 
lAlovetalJ i2004) for recent relativistic hydrodynamic simulations 
of such a scenario). Correspondingly, in the "collapsar" scenario, 



the core of a massive rota ting star collapses to form a rotating 
black hole lWoosle^l200lh which subsequently grows accreting 
the surrounding envelope of the star through a geometrically thick 
disc. An energetic, long-duration (i.e. > 2 s) GRB accompanied 
by a Type Ib/Ic supernova may be produced as a result of such 
(thick) disc accretion process. As in the neutron star merging sce- 
nario, part of the large amounts of the energy released by accretion 
is deposited in the low-density region along the rotation axis of 
the star, where the hea t ed gas expands in a jet-li ke fireball (see 
iMacFadven & Wooslevl 1 19991) : lAlov et alj feOOOl) for numerical 
simulations of this process). 

Even in the absence of other potential instabilities triggered 
by nonaxisymmetric perturbations or magnetic fields, the accretion 
torus that forms in the aforementioned catastrophic events is likely 
to be accreting on to the black hole with rates which largely ex- 
ceed the Eddington rate, a property which may lead to the rapid 
destruction of the disc. S uch a possibility was first suggested by 
lAbramowicz et all ll983l) . who pointed out how the equilibrium 
solution of a barotropic fluid orbiting a black hole with a general, 
non-Keplerian rotation law might be hydrodynamically unstable. 
The outermost closed equipotential surface of the torus has a dis- 
tinctive cusp at the inner edge through which mass transfer can be 
driven by small changes in the pressure gradient there (in a sta- 
tionary situation the pressure gradient at the cusp would be exactly 
zero). The instability is driven by the progressive penetration of the 
cusp into the disc as a result of the accretion of mass and angular 
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momentum and of the corresponding increase of mass and spin of 
the black hole. During the development of the instability the mass 
flux increases exponentially and leads to the complete disappear- 
ance of the disc into the black hole on a dynamical timescale, jus- 
tifying why it is usually referred to as the "runaway instability". 

Axisymmetric, general relativistic hydrodynamical simu- 
lation s of this process h ave been pe r forme d only re cently 
iFont & Daignd l2002albt IZanotti et ail 120031: iDaigne & Fontl 
2004). In these works, the disc is assumed to be non-selfgravitating 
and the relativistic hydrodynamics equations are then evolved in 
time in the curved spacetime of the black hole. The latter is in- 
stantaneously modified to new stationary solutions as the mass and 
spin of the black hole gradually increase as a result of accretion. 
While this is clearly an approximation to solving the full Einstein 
equations, it is a reasonable one when the disc-to-hole mass ratio is 
small enough and the changes in the spacetime between two succes- 
sive timelevels is correspondingly small. Within this approximation 
the numerical simulations have shown that while thick discs with 
a constant distribution of the specific angular momentum are sub- 
ject to the runaway instability irrespective of the initial conditions 
iFont & Daigner2 002a ; Zano tti et al . 2003), a small departure from 
a constant distribution of angular momentum (namely, a power- 
law increase with the radial distance with a small index) suffices to 
prevent the development of the instability iFont & Daignel f2002b ; 
Daigne & Font 2004). This result is in agreement with perturbative 
calculations of radially oscillating tori which indicate that the os- 
cillations have eigenfunctions that are significantly modified at the 
inner edge when the s pecific angular momentum is not constant 
iRezzolla et"ai]l2003bl) . However, no definitive conclusion can be 
made yet on the occurrence of the runaway instability as there are at 
least two other physical processes that could play some role. One of 
these is provided by the self-gravity of the disc which, from studies 
based on sequ ences of stationary models has been shown to favour 
the instability iMasuda et aill 998): the other one is brought up by 
the presence of magnetic fields that would affect globally the prop- 
erties of accretion flows and the angular momentum transport in 
discs surrounding black holes. 

Magnetic fields are indeed expected to influence significantly 
the physics of accretion discs, either geometrically thin or thick, in 
many respects, well beyond the issue of the runaway instability. In 
particular, the impact of magnetic fields in the dynamics of accre- 
tion flows has been recognized in a large number of works, but only 
recently the first magnetohydrodynamical simulations i n the r ela- 
tivistic regime have become available (but see Yokosawa 1 1995) for 
earlier investigations). A number of authors I McKinnev & Gammie 
| 2002t iDe Villiers et aljf2003t iGammie. McKinnev & TotrJ 120031 
Iffirose et alj 120031: iGammie. Shapiro & McKinnev! 120041) have 
started programs to perform realistic simulations of magnetically 
driven accretion flows in rotating black hole spacetimes. A funda- 
mental idea under intense scrutiny is that the combined presence of 
a magnetic field and of a differentially rotating fluid, in particular 
a Keplerian disc, can le ad to the so called magnetorotational insta- 
bility (see Balbus 1 2003)) for a review), generating an effective vis- 
cosity and thus transporting angular momentum outward through 
magnetohydrodynamic turbulence. 

Notwithstanding that magnetic fields play an important role 
in accretion theory, two new interesting results have emerged re- 
cently in purely hydrodynamical axisymmetric simulations. Both 
of them refer to the evidence that upon the introduction of per- 
turbations stable relativistic tori manifest a regular, long-term os- 
cillating behaviour. This, in turn, is responsible for the periodic 
changes of the mass quadrupole moment of the discs, thus mak- 



ing these objects promising new sources of detectable gravitational 
radiation. This is particularly true when the avera ge density in the 
discs is close to nuclear matter density IZanotti et alJl2003f) . as it 
is the case in discs formed as a result of binary neutron star co- 
alescence. In addition, when the thick discs are composed of the 
low-density material stripped from the secondary star in low-mass 
X-ray binaries (LMXBs), their oscillations could help explaining 
the high frequency quasi-periodic oscillations (HFQPOs) observed 
in the spectra of X-ray binaries iRezzolla et al. 2003b). In a model 
proposed recently and based on the evidence that tori around black 
holes have the fundamental mode of oscillation and the first over- 
tones in the harmonic sequence 2 : 3 : 4 ... to a good precision 
and in a very wide parameter space, the HFQPOs are explained in 
terms of p-mode oscillations of a small-size accretion torus orbiting 
around the black hole (see also Abramowicz & Kluzniak 2004 for a 
recent review on the phenomenology associated with HFQPOs and 
for alternative models). 

The aim of this paper is to extend the analysi s of the dy- 
namic s of axisymmetric relativistic tori carried out by Zano tti et alJ 

(2003) (henceforth paper I), in which the attention was limited to 
models with constant specific angular momentum orbiting around 
a Schwarzschild (nonrotating) black hole. Now we consider the 
more general case of nonconstant angular momentum thick discs 
in the rotating background spacetime provided by the Kerr m etric. 
Furthe rmore, since the results reported recently bv lDaigne"& Font 

(2004) pose limits on the occurrence of the runaway instability in 
the case of non self-gravitating, nonconstant angular momentum 
discs, we do not investigate this issue further here. Instead, we fo- 
cus on the study of the oscillating behaviour of marginally stable 
discs. To this aim we combine two complementary tools: nonlin- 
ear hydrodynamic simulations and a linear perturbative approach. 
With the first one we can follow the long-term evolution of per- 
turbed thick accretion discs, analyze their dynamics, and compute 
the gravitational wave emission in the Newtonian quadrupole ap- 
proximation. Correspondingly, with the second approach we can 
investigate how the oscillation properties of these objects depend 
on their geometrical features, surveying a very large parameter 
space. When applied to vertically integrated relativistic tori, this 
second approach has already shown its utility for explaining the re- 
su lts of the numerical simulations of Font & Daigne 1 2002a b|) and 
of lZanotti eta?]<2003l) . 

The paper is organized as follows: In Section|2|we briefly re- 
view the main properties of relativistic tori which obey a power- 
law angular momentum distribution. Next, in Section|3| we present 
the two mathematical approaches used in our investigation, namely 
the techniques for the numerical solution of the hydrodynamic 
equations and the basic aspects related with the solution of the 
eigenvalue problem in the linear perturbative framework. Section|4] 
presents the initial models, while Section [5] contains the discus- 
sion on the results. Finally, Section|6|is devoted to the conclusions. 
Throughout the paper we use a space-like signature (— , +, +, +) 
and a system of geometrized units in which G = c — 1. The unit 
of length is chosen to be the gravitational radius of the black hole, 
r g = GM/c 2 , where M is the mass of the black hole. Greek inde- 
ces run from to 3 and Latin indeces from 1 to 3. 



2 STATIONARY FLUID CONFIGURATIONS 

The procedure for building a stationary and axisymmetric fluid con- 
figuration orbiting around a Kerr black hole and obeying a power- 
law distribution of the specific angular momentum in the equatorial 
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plane has been reviewed in great detail by Daig ne & FonJ 120041) . 
Here, for the sake of completeness, we will only recall those defini- 
tions and physical properties that are essential to the present study. 
In what follows we consider a perfect fluid described by the stress- 
energy tensor 

T' 1 " = (e + p)u Al M" + pg^ v — phu^u" -\-pg^ v , (1) 

where g l>v are the coefficients of the Kerr metric in Boyer- 
Lindquist coordinates (t, r, 9, (j>) and e, p, p, and h = (e + p)/p 
are the energy density, the isotropic pressure, the rest-mass density, 
and the specific enthalpy, respectively, each of them measured in 
the frame comoving with the fluid. The fluid is supposed to obey a 
polytropic equation of state p = up 1 (EOS), where k is the poly- 
tropic constant and 7 is the adiabatic index. We use fl = u 9 /u 1 to 
denote the angular velocity of the fluid and I = —u^/ut to denote 
its specific angular momentum. Moreover, the rotation law on the 
equatorial plane is given by a power-law distribution of the specific 
angular momentum 

£(r, 6 = tt/2) = Sr q , (2) 

where 5 can be either positive or negative, according to the disc ro- 
tation being prograde or retrograde, respectively, with respect to the 
black hole rotation. When the motion of the fluid is just circular and 
does not change in time (i.e. dt = = d$), the relativistic Euler 
equations in the r and 6 directions are given by the (Bernoulli-type) 
equilibrium conditions 

-_ = - Vi ln( Wt ) + T -^ (, = r,0), (3) 

where the right-hand-side is nothing but the opposite of the rela- 
tivistic four acceleration a% = u a \7 a Ui, which vanishes in the case 
of a purely geodetic motion 1 . The procedure for solving Eq. ^3} ex- 
ploits the fact that the surfaces of constant Q, (the so-called von 
Zeipe l cylinders) coi ncide with the surfaces of constant £ (see 
iDaigne & Font! l2004[) for details), a res ult which holds true only 
for a barotropic fluid fAbramowi czll97ll) . One of the most relevant 
features of the equilibrium solution, irrespective of the distribution 
of the angular momentum, is that the resulting thick disc can have 
two Keplerian points on the equatorial plane (i.e. points where the 
rotation law is that of the Keplerian motion): the first one is the 
cusp, through which matter can accrete onto the black hole, and 
the second one marks the position of the maximum rest-mass den- 
sity. However, while in the case of a constant angular momentum £ 
the position of the cusp r cusp is always smaller than the marginally 
stable orbit r ms , when the angular momentum is not constant, it 
may also happen that r cusp > r ms , depending on the choice of the 
parameters 5 and q in Eq. J2j. 

Among the large family of initial configurations that are in 
principle possible, we only consider those that possess a cusp, a 
centre, and with a closed equipotential through the cusp. This only 
happens when < q < 0.5 and when |<S m8 | < |<S| < |<SmbL where 
Sms a nd 5 m b are two limiting values defined in Daigne & Font 
(2004) (the subindex 'mb' refers to the marginally bound orbit). 
Also note that relativistic tori with £ following a power-law gener- 
ally have larger sizes than constant angular momentum tori, a fact 
which can produce initial models with central densities one or two 
orders of magnitude smaller. 

1 Note that Eq. |3j implies that a stationary and axisymmetric circular mo- 
tion of a perfect fluid is geodesic if and only if it follows a Keplerian rotation 
law. 



3 MATHEMATICAL FRAMEWORK 

The results we report in this paper are based on two complemen- 
tary approaches involving both nonlinear hydrodynamic simula- 
tions and the solution of the eigenvalue problem within a pertur- 
bative linear analysis. We recall in what follows the basic aspects 
related to both of them. 



3.1 Solution of the nonlinear hydrodynamics equations 

The nonlinear, axisymmetric, general relativistic code used in the 
si mulations we report is an extended version of the one presented 
bv lZanotti et al J 120031 which incorporates the effects of a rotating 
background spacetime. In this code the hydrodynamic equations 
are implemented as a first-order, flu x-conservative system accord- 
ing to the formulation developed by Banvul s"et alJ 1199% . The ex- 
plicit form of thi s system of eq uations for the Kerr metric can be 
found in Font & Daigne 1 2002a). The equations are solved using a 
high-resolution shock-capturing scheme based on an approximate 
Riemann solver (either HLLE or Marquina's). Second-order accu- 
racy in both space and time is achieved by adopting a piecewise- 
linear cell reconstruction procedure and a second order, conserva- 
tive Runge-Kutta scheme, respectively. 

The computational grid consists of 7V r x N$ zones in the ra- 
dial and angular directions, respectively. The innermost zone of the 
radial grid is placed at r m i n = rhorizon + e, while the outermost 
radial grid point is at a distance about 30% larger than the outer 
radius of the torus, r out . While e depends on the particular model 
considered, it is typically 0.1r fl or smaller. The radial grid is built 
by joining smoothly two patches obtained using two different algo- 
rithms. The first part extends from r m i n to the outer radius of the 
torus. It is logarithmically spaced and has the maximum radial res- 
olution at the innermost grid zone, Ar = 1 x 10~ 3 . The spacing 
of the second part of the radial grid is uniform and it extends up to 
r mal . The total number of radial grid points N r is selected so as to 
have a radial resolution at the outer boundary Ar ~ 0.2. Typically, 
we use iV r ~ 250 for the more compact tori (radial length L < 15), 
while N r ~ 350 for the more extended models. The angular grid, 
on the other hand, which consists of Ne = 84 zones in all of the 
simulations, covers the domain from to n. 

As customary for finite difference hydrodynamical codes 
which cannot handle vacuum regions, a low density "atmosphere" 
is introduced in those parts of the numerical domain not occupied 
by the torus. In all our simulations we have chosen a special versio n 
of the spherically accreting solution described by Mi chell ll972h . 
suitably modified to account for the rotation of the black hole (see 
Appendix A for a detailed description). Since this atmosphere is 
evolved as the rest of the fluid, one should simply take care that its 
dynamics does not interfere with that of the object being studied. 
We have verified that this is indeed the case if the maximum den- 
sity of the atmosphere is 5 — 6 orders of magnitude smaller than 
the central density of the torus. 

As mentioned in the Introduction we are not interested here 
in the study of the runaway instability but rather we focus on the 
oscillation properties of thick discs filling their outermost closed 
equipotential surface (we also refer to these discs as "marginally 
stable"). Hence, we assume that the background spacetime is sim- 
ply the one provided by the Kerr metric and that it does not change 
during the evolution; this pr events the dev elopment of the insta- 
bility iFont & Daignell2002j:IZanotti et alJ2 003) by construction. 
This is a reasonable approximation in the present context for two 
different reasons. The first one is that most of the tori we evolve 
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have rest masses M t which are much smaller than that of the black 
hole (i.e. Mt/M — 0.1). Hence, we can neglect their contribution 
to the overall gravitational field when compared to that created by 
the black hole. The second reason is that a small power-law index 
q in the angular momentum distribution is enough to reduce sig- 
nificantly the amount of rest-mass accreted onto the black hole so 
that, effectively, the mass and spin of the black hole do not change 
significantly over the timescale of our s imulations . This was shown 
numerically by Font & Daigne 1 2002b), Daiene & Font 1 2004) and 
explained in the perturbative analysis by |R£zzcJki^HdTl20o"3bl) . To 
validate this approximation we monitor the rest-mass accreted by 
the black hole according to the formula 

m(r min ) = -2tt / ^ P Wv r d9\ , (4) 

JO lr min 

where g is the determinant of the metric and where r m i n is the ra- 
dius of the innermost radial cell. Similarly, we monitor the angular 
momentum flux across the horizon as 

j(r mi „) ee -2tt / l^g~pWv r de\ . (5) 

Note that W and v r in Eqs. J4j and J5J are the Lorentz factor and 
the coordinate radial velocity, respectively, as measured by the Zero 
Angular Momentum Observer. In all of the simulations performed 
here the global relative changes in the black hole mass and spin 
resulting from accretion are AM/M ~ AJ/J < 10~ 4 . Hence, 
the assumption of a fixed spacetime is accurate and justified. 

3.2 Solution of the eigenvalue problem 

In addition to the hydrodynamic simulations we also perform a lin- 
ear perturbative analysis of the axisymmetric modes of oscillation 
of relativistic tori in a Kerr spaceti me. This is done following the 
procedure outlined in iRezzolla et ail §003b ) for the Schwarzschild 
spacetime and recently extended to the Kerr case by Montero et al. 
(2004). The method, whose details can be found in the aforemen- 
tioned references, consists in solving the perturbed relativistic con- 
tinuity and Euler equations obtained after introducing perturbations 
with a harmonic time dependence in the velocity and the pressure. 
The system of perturbed hydrodynamical equations is then cast 
into a set of coupled ordinary differential equations and solved as 
an eigenvalue problem, where the perturbed quantities are treated 
as eigenfunctions and where the eigenvalues provide the eigenfre- 
quencies of the s ystem. A similar p rocedure but for thin discs has 
been described by Rodriguez et al. (2002). 

Two approximations are adopted in the perturbative analysis. 
First, we neglect the perturbations of the background spacetime, 
something which is usually referred to as the Cowling approxima- 
tion lCowling|l94ll) . This choice is consistent with that of neglect- 
ing the contribution of mass and angular momentum fluxes to the 
mass and spin of the black hole in the nonlinear hydrodynamic sim- 
ulations. The second approximation is that of adopting a vertically 
integrated description of the tori. While this approach simplifies 
the numerical treatment considerably (the angular dependence is 
removed and only the radial one remains) and allows for a satisfac- 
tory calculation of the eigenfunctions and eigenfrequencies for the 
lowest-order modes, it may not be sufficiently accurate to capture 
the properties of the higher-order modes, that intrinsically have a 
larger number of nodes in the radial and polar directions. This ap- 
proximation produces some discrepancies between the two comple- 
mentary approaches we use, for instance, when exciting selectively 
higher-order modes, or when comparing the eigenfrequencies of 
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Figure 1. Distribution of the whole set of equilibrium models in the plane 
(fcusp — fms, L). Notice that models Alb and Ale (not plotted) share the 
same location of model Ala. 

modes larger than the first overtone. A more detailed discussion of 
this issue is presented in Section l5.1.2l . 



4 INITIAL MODELS 

The initial models for the numerical simulations are a sequence 
of relativistic tori which fill their outermost closed equipotential 
surface and thus their inner radius coincide with the position of 
the cusp, i.e. r- m = r cusp . TableQreports the main properties of 
the unperturbed models. All of the models but two (i.e. A5b and 
A5c) are built with an adiabatic index 7 = 4/3. More specifically, 
models A5b and A5c have 7 = 5/3 and 7 = 2, respectively. 

The value of the constant 5 in the power-law distribution, 
Eq. 0J> spans the allowed range between <S ms and 5 m b, and is se- 
lected via the additional parameter A as S — min(S m s, <S m b) + 
A|5 m b — S ms |- The resulting value of 5, together with the power- 
law index q, are the most important factors in determining the geo- 
metrical properties of the torus, in particular for the location of the 
Keplerian points (i.e. r cusp and r max , not reported in TableQ and 
for the radial size of the disc L = r out — ri n . 

Figure Q provides a graphical view of the distribution of the 
initial models in the plane (r cusp — r ms , L). We recall that the po- 
sition of the marginally stable orbit, r ms , is a decreasing function of 
the black hole spin parameter a, while by increasing the power law 
index q the inner radius of the tori gets larger. As a result it is a natu- 
ral consequence of nonconstant angular momentum tori in the Kerr 
metric to have r cusp — r ms usually larger than zero, making accre- 
tion more difficult than in the case of constant angular momentum 
tori. Only for small values of a and q can the torus penetrate deeper 
in the potential well, pushing the cusp below the marginally stable 
orbit. 

It is worth underlining that, once perturbed, the numerical 
evolution of the stationary models reported in Table Q can indeed 
provide an insight into the dynamics of astrophysical thick discs 



Oscillating relativistic tori around Kerr black holes 5 



Table 1. Fundamental properties of the initial models. From left to right the columns report the name of the model, the black hole spin parameter a, the 
torus-to-hole mass ratio Mt / M, the power-law index q of the angular momentum distribution, the parameter A fixing the value of the constant <S, the adiabatic 
index 7, the polytropic constant k, the inner and the outer radius of the tours, r; n and r out , the equatorial size of the torus L, the orbital period at the point of 
maximum rest-mass density t or t, , the maximum rest-mass density, and the average rest-mass density. For all models considered the mass of the black hole is 

M = 2.5M . 
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3.19X10 13 


5.69X10 12 


D4a 
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0.1 


0.2 


4/3 


1.86xl0 13 


3.230 


7.202 


3.97 


0.91 


1.17X10 14 


3.37X10 13 


Ela 
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4/3 


4.12X10 14 


1.971 


16.535 


14.56 


0.71 


1.80xl0 13 


1.04xl0 12 
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4.65 xlO 14 


1.938 


25.801 


23.86 


0.75 
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2.74X10 11 


E3a 


0.9 


0.1 


0.2 


0.3 


4/3 


3.62 xlO 13 


2.611 


8.270 


5.66 


0.85 


9.38X10 13 


1.69xl0 13 


E4a 


0.9 


0.1 


0.2 


0.4 


4/3 


7.16xl0 13 


2.519 


10.459 


7.94 


0.91 


4.58 xlO 13 


6.74 xlO 12 


E5a 


0.9 


0.1 


0.2 


0.5 


4/3 


1.26xl0 14 


2.443 


13.500 


11.05 


0.98 


2.37X10 13 


2.69X10 12 


A5b 


0.0 


0.1 


0.2 


0.2 


5/3 


1.28xl0 9 


6.824 


17.118 


10.29 


2.97 


9.25 xlO 12 


3.22X10 12 


A5c 


0.0 


0.1 


0.2 


0.2 


2 


9.82xl0 4 


6.824 


17.118 


10.29 


2.97 


7.18X10 12 


3.21 xlO 12 



*Note that for model Ala the constant A is not defined and I = 3.8. 



around stellar-mass black holes. All of the astrophysical scenar- 
ios where relativistic tori are expected to form (such as the merger 
of two neutron stars, or of a black hole and a neutron star, or in 
the core collapse of a massive star) clearly suggest that the newly 
formed torus will be far from equilibrium and may indeed represent 
only a transient phase to a process which will ultimately destroy the 
torus. 

For this reason we consider three different types of initial per- 
turbations, aimed at reproducing some of the physical conditions 
that might occur in realistic astrophysical scenarios. The first one 
has already been described and motivated in paper I and is based 
on the inclusion of a radial velocity (we recall that in equilibrium 
all the velocity components but the azimuthal one are zero). This 
perturbation is parametrized in terms of a dimensionless coefficient 
77 of the sphericall y symm etric accretion flow onto a Schwarzschild 
spacetime lMichei ri972l) . i.e. v r = ?7(ur)Michci- The second type 
of perturbation intends to affect the rest-mass density only, either in 
a global way (e.g. by rescaling the stationary profile), or in a local 
random way (e.g. by introducing small variations randomly dis- 
tributed within the torus). Finally, the third type of perturbation is 
closely related to the perturbative analysis discussed in Section lT2l 
and uses the eigenfunctions of the p-modes obtained in the linear 
analysis as initial data for the perturbations in the density and veloc- 
ity. This procedure is particularly effective when the lowest-order 
modes are used. 



5 NUMERICAL RESULTS 
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Figure 2. Time evolution of the mass accretion rate for representative tori 
orbiting around a Schwarzschild black hole. The figure shows the stabilizing 
effect of increasing the power law index from q = (model Ala) to q = 
0.1 (model A3a), which reflects in progressively smaller accretion rates. In 
all of the models the perturbation factor r\ = 0.08. 



The disc-to-hole mass ratio considered in our initial models 
(M t /M = 0.1) as well as the power-law index of the angular 
momentum distribution chosen, make these objects stable with re- 
spect to the runaway instability. Indeed, as reported recently by 



iDaiene & Fontl 120041) for discs with M t /M = 0.1, the critical 
power-law index separating stable and unstable models appears to 
be q CT ~ 0.05 — 0.06 for mass accretion rates of about a few times 
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MqsT 1 . Hence, when perturbed, these models will only respond 
with oscillations that are harmonic as long as the amplitude of the 
perturbation is sufficiently small (see paper I for a discussion of the 
transition to nonlinear oscillations). As a result, with the exception 
of model Ala, which has a power law index q = 0, we expect the 
models selected in Table Q to be particularly suitable for studying 
their response to perturbations and, in particular, to compute the 
associated mode eigenfrequencies (Section l5.1> and gravitational 
wave emission (Section f5.2t . 

5.1 Oscillation properties 

5.1.1 Dynamics of representative models 

Figure[2]shows the time evolution of the accretion rate of the first 
three models of TableQ namely models Ala, A2a, A3a, all of them 
referring to a Schwarzschild black hole. The time evolution dis- 
played in Fig.|2|corresponds to models having an initial perturba- 
tion in the radial velocity with an amplitude r\ = 0.08. The final 
time plotted in the figure corresponds to 100 orbital periods, which 
are measured with respect to the orbital period of the rest-mass den- 
sity maximum in the unperturbed torus. 

From model Ala to model A3a the power law index changes 
from to 0.1, thus providing a progressive deviation from a con- 
stant angular momentum distribution toward a steeper power law 
distribution. As it is obvious from the figure when the power law 
index q increases while maintaining A constant, the mass accretion 
rate decreases correspondingly. Already at q = 0.05 (model A2a) 
the behaviour of rn is quite different from that of a constant angular 
momentum model, and at q = 0.1 the behaviour of the mass flux is 
entirely dominated by the behaviour of the accreting atmosphere. 

We have analyzed the time evolution of the accretion rate for 
the rest of models reported in Table Q finding that apart from an 
early intense burst of accretion due to the response of the system to 
the initial perturbation (this is evident e.g. in model A3a in Fig|2j, 
none of them shows a significant variation of the position of the 
inner radius when compared to the initial one, thus having an es- 
sentially suppressed mass flux. This behaviour does not depend in a 
substantial manner on whether the inner radius is smaller or larger 
than the marginally stable orbit. This is the case, for instance, of 
model Cla which although has r cusp — r ms ~ —1.0 is very stable 
and barely accretes on to the black hole, with m ~ 10 _6 Mqs _1 . 
On the contrary, this behaviour strongly depends on the value of the 
power law index q. This reflects the fact that tori with nonconstant 
distributions of specific angular momentum have oscillations which 
are significantly damped at the inner edge (Rezz ollaetall2003bh . 
We also recall that the amplitude of the accretion rates reported in 
Fig.[2]depends on the perturbation factor 77, and that, as shown in 
paper I, this dependence is linear for 77 < 0.04 (cf. Fig. 9 of paper 
I). 

Despite a distinctive quasi-periodic behaviour is quite appar- 
ent from Fig.|2| we note that for most models the evolution of the 
accretion rate reflects mainly the response of the atmosphere to the 
perturbations propagating from the torus. For this reason the peri- 
odic character of the dynamics of the discs becomes more evident 
in the time evolution of the maximum rest-mass density rather than 
in the mass flux. This is reported in Fig. [3] for four representative 
models, two of them corresponding to a Kerr black hole of spin 
rate a = 0.7 (Dla, D2a), and the other two corresponding to a Kerr 
black hole of spin rate a — 0.9 (E3a, E4a). After the tori relax 
from the initial perturbation at 7j/i or b ~ 2, they start oscillating at 
regular, quasi-periodic intervals. This is a feature which discs with 



Table 2. Frequency ratio o\ // as extracted from the power spectra of the 
L/2 norm of the rest mass density (second column) and as computed from 
the solution of the eigenvalue problem (third column). 



Model (oi/f) num (oi/f)i i: 



Ala 


1.45 


1.47 


A2a 


1.46 


1.41 


A4a 


1.36 


1.36 


A5a 


1.31 


1.26 


Dla 


1.37 


1.32 


D2a 


1.40 


1.36 


D3a 


1.35 


1.28 


D4a 


1.35 


1.28 


Ela 


1.33 


1.40 


E3a 


1.31 


1.25 


E4a 


1.36 


1.29 


E5a 


1.38 


1.32 



power-law distributions of angular momentum share with constant 
angular momentum discs as those investigated in paper I. The har- 
monic variation of the hydrodynamical quantities is analyzed in 
more detail next. 

5.1.2 Comparisons with the perturbative analysis 

Additional information on the quasi-periodic behaviour of the hy- 
drodynamics variables discussed in the preceding section can be 
extracted through a Fourier analysis of the corresponding time evo- 
lutions. For this purpose we have calculated the Fourier transforms 
of the time evolution of the L2 norm of the rest-mass density for 
all models, defined as | \p\ | 2 = E^fi Pij- As this is a g lobal 

quantity it is particularly useful for comparisons with the results 
of the eigenvalue problem in the linear perturbative analysis. Fur- 
thermore, we can also compare our findings with those reported in 
paper I for the case of constant angular momentum discs. 

A first result that emerges from the Fourier analysis is that the 
overall dynamics of nonconstant angular momentum discs is more 
complex than in the case of constant distributions of the angular 
momentum. This is reflected in particular in the fact that the power 
spectra of the L2 norm of the rest-mass density show, in general, 
a richer structure. This can be seen in Fig. [4] where we plot the 
power spectra of the L2 norm of p, obtained with a Hanning filter 
(Press et al. 1996) for one representative model, namely E3a. 

Although not all of the computed spectra show the features of 
Fig. [4] with the same clarity, all of them present the same essen- 
tial characteristics, namely a fundamental mode / and a series of 
overtones. Among them we can recognize the two first overtones 
predicted by the linear analysis and indicated as o\ and 02 . In ad- 
dition to these, further modes appear as linear combinations of /, 
01 and 02, the most prominent of which is the one at 2/. Since 
the ratio o\ // can be computed with fairly good accuracy from the 
spectra, we compare it with the ratio between the two first eigenfre- 
quencies computed with the linear perturbative analysis. This com- 
parison is reported in Table|2| Note that, for those models for which 
it was possible to identify the mode o\ unambiguously, the agree- 
ment with the prediction of the linear perturbative approach is very 
good, with differences of ~ 5% at most. A similar comparison 
can also be done for the mode 02, although its identification from 
the computed spectra is much more difficult. In the two cases when 
this was possible, namely models E3a and D4a, the ratio 02/ f was 
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Figure 3. Time evolution of the central rest-mass density normalized to its initial value for models Dla, D2a (black hole spin a = 0.7) and E3a, E4a (black 
hole spin a = 0.9). The perturbation factor j) = 0.08 for models of both panels. 
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Figure 4. Power spectrum of the L2 norm of the rest-mass density for 
model E3a. The units in the vertical axis are arbitrary and the power spec- 
trum was obtained using a Hanning filter. 



shown to differ from the one predicted by the linear analysis by less 
than 7%. 

It is worth emphasizing that the excitation of a particular over- 
tone depends sensibly on the choice of the initial perturbation. In 
particular, we can excite selectively a specific mode by perturbing a 
given equilibrium model through the use of the vertically integrated 
eigenfunctions for the velocity and rest-mass density that have been 
calculated for that model through the perturbative analysis. Exploit- 
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Figure 5. Power spectra of the L2 norm of the rest-mass density for model 
E3a. The different lines refer to different initial perturbations: the solid line 
corresponds to an initial perturbation given by a parametrized spherically 
symmetric radial velocity, the dotted line to a global initial perturbation 
in the density, and the dashed line to an initial perturbation given by the 
eigenfunction of mode 01 . Vertical units are arbitrary. The values have been 
rescaled in order to match the power of the fundamental frequency. 



ing this possibility we have investigated in detail the dynamics of 
model E3a, considering initial data which consist either of generic 
perturbations in the radial velocity (with r\ — 0.08) or in the rest- 
mass density (or in both), or of specific perturbations obeying the 
eigenfunctions for the velocity and the rest-mass density for the 
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oi overtone (the fundamental mode also would be excited in this 
case). The results of this investigation are summarized in Fig. [5] 
which shows the power spectra of the L2 norm of p for the three 
different initial perturbations. As it is clear from this figure the two 
generic initial perturbations, represented by a global perturbation 
in the radial velocity (solid line) and by a global perturbation in 
the rest-mass density (dotted line), are less efficient in exciting the 
mode 01 than when the mode eigenfunction is used (dashed line). 
In the latter case the power channeled in that mode increases by a 
factor ~ 20 with respect to initial perturbations in the global veloc- 
ity and by a factor ~ 5 with respect to the initial perturbations in 
the global density. 

We note that exciting overtones of the fundamental mode 
above o\ is increasingly difficult. The tests that we have performed 
using the eigenfunctions of the mode 02, for instance, could not 
provide a clear signature of a selective excitation, in contrast with 
the mode 01 . We believe this is probably due to the approximation 
made in the linear perturbative approach in treating thick discs as 
vertically integrated objects. While we expect this approximation to 
be a satisfactory one for the lower-order modes, it becomes gradu- 
ally less accurate as the mode number increases and the small-scale 
features of the eigenfunctions become more important. 

We have not yet commented on the presence of the modes as 
linear combinations of /, 01 and 02 shown in Fig.[4] and in partic- 
ular on the overtone appearing at twice the fundamental frequency, 
which is, after the fundamental mode, where most of the power is 
concentrated. A plausible interpretation of the peak at 2/, and of 
the other linear combinations shown in the spectrum of Fig. [4] is 
that they are the result of a nonlinear coupling effect. It is, in fact, 
a general property of nonlinear systems in the limit of small oscil- 
lations that of showing nonharmonic oscillations, i.e. linear com- 
binations of their normal modes of oscillations. Thus, if the sys- 
tem has eigenfrequencies Ui, the nonlinearity of the equations will 
also produce modes at frequencies uji ± uij (cf. Landau & Lifschitz 
1976, §28), with amplitudes which are proportional to the product 
of the amplitudes of the combining frequencies. It should be noted 
that this nonlinear coupling is particular evident in the present cal- 
culations which have been performed with an initial perturbation 
amplitude 77 = 0.08, at which the system's response is already 
nonlinear 2 . 

A consequence of this nonlinear coupling among modes is 
that, if Af and A2] denote the power spectra amplitudes of a 
generic quantity A (e.g. the rest-mass density) at the frequencies 
/ and 2/, respectively, we then expect that Af S> Aif and the 
two amplitudes to scale like Aif oc Aj. We have tested that this 
is indeed the case by considering one representative model, D2a. 
For this model we have computed from the power spectrum a se- 
quence of pairs (log A / , log Aif ) for different values of the per- 
turbation parameter r\, verifying that to a good approximation they 
settle along a strai ght line of s lope 2 in the corresponding plane. 

In|Rezzolla etai] fe003bh it was shown through linear analysis 
that p-modes in a thick disc with a constant distribution of angular 
momentum have frequencies which stay in the harmonic sequence 
2:3:4.... Since this same sequence was also found in the non- 
linear numerical simulations reported in paper I it was thus natural 
to interpret them as an evidence of the excitation of p-modes. In 
particular, the mode at twice the fundamental frequency was iden- 
tified as the mode 02 , for which 02/ f = 2. However, it is now clear 



The maximum perturbation amplitude triggering linear perturbations has 
been estimated to be 77 ~ 0.04 in paper I (cf. Fig. 9 of that paper). 
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Figure 6. Ratio of overtones as computed through the linear analysis. The 
ratios ol/f and o2// approach 3/2 and 2, respectively, as the cusp of the 
torus penetrates deeper in the potential well. Three of the models considered 
here (geometric symbols) belong to the sample of TablelTI For comparison, 
the filled geometric symbols report the values obtained from the numerical 
simulations. 



that in the case of constant angular momentum discs the peak at 2/ 
receives contributions both from the mode 02 and from the nonlin- 
ear coupling / + /. As a result, part (if not most) of the power at the 
peak at twice the fundamental frequency in Fig. 7 of paper I could 
be the result of the nonlinear coupling discussed here. 

More recently, iMontero et alj |2004) have investigated p- 
modes in thick discs with more generic angular momentum dis- 
tributions and found that in this case deviations from the simple 
sequence 2:3:4... can be large. We show in Fig. [6] the fre- 
quency ratios o%/f and 02/ f for a number of models of our sam- 
ple as a function of the penetration of the inner radius of the discs 
in the black hole potential well. All of the values in Fig. [6] have 
been computed with the linear code and some of the unperturbed 
models coincide with those evolved with the nonlinear code. These 
are E3a, E4a, E5a and are indicated with the geometric symbols. 
Note that both 01/ f and 02 // differ from 3/2 and 2, values which 
they approach only in the limit of a torus penetrating deeply in the 
potential well. Furthermore, while the spectra obtained from the 
numerical simulations show a clear peak at the position predicted 
by the linear analysis for the mode 01, it is difficult to find a peak at 
the position predicted by the perturbative code for the 02 overtone. 

In summary, the comparison between the results of the non- 
linear numerical simulations and those coming from the linear per- 
turbation analysis indicate that the fundamental mode / and its first 
overtone in the numerical simulations do represent the first two p- 
modes of the system and that these are in ratio 01// close to 3/2, 
with deviations that can be as large as ~ 15%. Additional p-modes 
are also probably excited in the simulations, namely the mode 02 
as shown in Fig.|4]for model E3a, but the corresponding power is 
in general too small to be visible in the spectra of the rest of the 
models. On the other hand, the computed spectra show overtones 
at integer multiples of the fundamental frequency (most notably at 
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2/), plus additional modes which are linear combinations of / and 
ox, and that are all the result of nonlinear coupling effects. 

The present results also s hed some addit i onal l ight on the 
physical mechanism proposed by Rezzolla et al. 1 2003a) to explain 
HFQPOs in black hole candidates. The main idea is that a partic- 
ular class of observations showing frequen cies in the ra tio 2 : 3 
and 1 : 2 iRemillard et all2002LrAbramowicz et all2003l) could be 
explained in terms of the excitation of the p-modes of a thick accre- 
tion disc orbiting around a black hole. According to the numerical 
simulations presented here, in the linear regime the 02 p-mode will 
be little excited, since very little power is in general channeled in 
this mode. However, a mode at 2/ will be present and strongly ex- 
cited in the nonlinear regime. As a result, the ratio 1 : 2 is satisfied 
with very good precision, is independent of the distribution of the 
angular momentum and appears as soon as the nonlinear coupling 
is triggered. The ratio of the fundamental mode and of the first over- 
tone, on the other hand, is more sensitive on the distribution of the 
angular momentum and may differ from 2 : 3 up to 15%. 

We conclude this section by recalling that the fundamental 
modes are closely related to the epicyclic oscillations at the lo- 
cation of the maximum rest ma ss densi ty of the disc and tha t, as 
shown by Rezzo lla et alJ l2003b) and by Montero et al. 1 200J), the 
eigenfrequency of the fundamental p-mode of the disc tends in the 
limit of vanishing size to the radial epicyclic frequency at the disc 
rest-mass density maximum. 



5.2 Gravitational wave emission 

The oscillating behaviour of the perturbed tori that we have dis- 
cussed in the previous section is responsible of significant changes 
of the quadrupole moment of these sources. When the tori are com- 
pact enough and dense enough (e.g. when they are produced via bi- 
nary neutron star mergers) the changes of the quadrupole moment 
result in the emission of potentially detectable gravitational radia- 
tion IZanotti et all2003l) . In paper I it was shown that the amount 
of energy that can be radiated in this way from perturbed "toroidal 
neutron stars" can indeed be very high, with resulting signal-to- 
noise (S/N) ratios comparable or larger than those expected from 
the burst signal associated with the bounce in the collapse of the 
core of massive stars (but see Miiller et al. 2004 for calculations 
of the gravitational radiation resulting from convection in the core 
collapse scenario). This property is closely related to the toroidal 
topology of these objects, which have their maximum rest mass 
density off-centered and thus have intrinsically high quadrupole 
moments and significant time variation of the latter. In the follow- 
ing we reexamine the issue of the gravitational wave emission from 
nonconstant angular momentum tori orbiting around Kerr black 
holes, and compare our findings with the results of paper I. 

The procedure for computing the gravitational wave emission 
is the same presented in paper I and it is based on the use of the 
Newtonian quadrupole approximation. In particular, the time evo- 
lution of the wave amplitude A20 , which is the second time deriva- 
tive of the mass quadrupole moment, is first computed through the 
"stress formula" iFinn 1989) to give 
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where k = 167T 3 / 2 /Vl5, z = cos 6 and $ is the gravitational 
potential, approximated to the first Post-Newtonian order as $ = 
(1 — <?rr)/2. The transverse traceless (TT) wave amplitude is then 
computed as (Zwerger & Miiller 1997) 



h TT (t) 



E2 



A 20 

R 



(7) 



where R is the distance to the source and where the detector's beam 
pattern function, F+ = F+ (R, 8, (/>), is set to one as if optimal con- 
ditions of detectability were met. In order to investigate the possi- 
bility that pulsating relativistic tori can be effectively detected by 
the interferometric instruments presently in operation or under con- 
struction, we have computed the signal-to-noise ratio with respect 
to LIGO I, Advanced LIGO, and VIRGO. For doing this we need 
an estimate of the frequency where most of the gravitational wave 
emission is channeled. This is provided by the characteristic fre- 
quency f c , which is a detector dependent quantity. If the detec- 
tor has a power spectral density Sh(f), then the characteristic fre- 
quency is given by <Thornell987t) 
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where (\h(f)\ 2 ), which is an average over randomly distributed 
angles of the Fourier transforms h(f) of h TT (t), has been approx- 
imated as (\h(f)\ 2 ) ~ \h{f)\ 2 - Strictly related to the character- 
istic frequency is the characteristic amplitude, that gives a typical 
measurement of the gravitational wave emission h TT detected by 
a given particular instrument and reads 
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The resulting signal-to-noise ratio is then computed as 
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where hims(fc) = \J fcSh(fc) is the strain noise of the detector 
at the characteristic frequency. 

FigureQ shows the strain noise of the three detectors consid- 
ered here as a function of frequency. Note that all of the sensitiv- 
ity curves reported assume an optimally incident wave in position 
and polarization (in agreement with the assumption of setting the 
detector's beam pattern function to unity) and that the curve for 
Advanced LIGO represents the sum of presently anticipated noise 
sources. Fig.|7|also shows the computed characteristic amplitudes 
in comparison with the strain noise of the detectors for sources 
located at a distance of 10 Kpc (i.e. for galactic sources) and of 
20 Mpc (i.e. for sources within the Virgo cluster). As it is evident 
from the figure, in the first case all of the models considered lie 
well above the sensitivity curves of the detectors, while in the sec- 
ond case only one model would be marginally detected by the Ad- 
vanced LIGO detector. Note also that the relatively wide scattering 
of the models in the figure (from the most powerful one, D4a, to the 
least powerful one, CI a) are due to the differences in the average 
rest-mass densities. 

Table|3|reports the characteristic frequency, characteristic am- 
plitude, and S/N ratio of the gravitational waves emitted by our 
tori models, with respect to the different detectors considered. In 
particular we choose the sources at the distance of lOKpc for LIGO 
I and VIRGO, while extragalactic sources at a distance of 20 Mpc 
for Advanced LIGO. Also for the case of nonconstant angular mo- 
mentum discs orbiting around Kerr black holes considered in the 
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Figure 7. Characteristic wave amplitudes for the tori models of Table Q 
with respect to the strain noise of LIGO I and the planned sensitivity of Ad- 
vanced LIGO, respectively. The amplitudes are computed at both a galactic 
distance of lOKpc and at an extragalactic distance of 20Mpc. The planned 
strain noise of VIRGO is also reported for comparison. 



present paper, Tableland Fig.[7]confirm the results found in pa- 
per I for constant angular momentum discs around Schwarzschild 
black holes, namely that oscillating relativistic tori can be promis- 
ing sources of gravitational wave emission. It is important to em- 
phasize that there are firm reasons to believe that the numbers 
reported in Table [3] should be considered as lower limits of the 
amplitudes of the gravitational wave signals that could be effec- 
tively emitted. Firstly, because the torus-to-hole mass ratio could 
be higher than 0.1, and the gravitational signal, which scales lin- 
early with this ratio (see paper I), would then be proportionally 
larger. Secondly, because the perturbations affecting the torus could 
be dramatically larger in a realistic formation scenario, i.e. after 
the core-collapse of a massive star or binary neutron star merger. 
Thirdly, because of issues having to do with the EOS of the material 
of the disc. We have actually verified that by changing the adiabatic 
index 7 of the polytropic EOS we use the gravitational radiation 
emitted can significantly increase. In particular, we have compared 
the wave amplitudes for models A5a, A5b, and A5c, which have 
7 = 4/3,5/3, and 2, respectively. While in the transition from 
7 = 4/3 to 7 = 5/3 the gravitational signal is essentially unmod- 
ified (it increases by ~ 3%), when the polytropic index is 7 = 2 
(stiff EOS) the gravitational signal increases by ~ 30%. Finally, 
we should also mention that Advanced LIGO has still plenty of 
improvement margins achievable which are currently under inves- 
tigation [for instance, by changing the position and the transmission 
of the signal recycling mirror I Shoemaker 2004)]. This provides an 
additional argument for considering the values reported in Table [5] 
as lower limits of the amplitudes of the gravitational wave signal 
from relativistic thick discs. 



6 CONCLUSIONS 

We have studied the dynamics and the gravitational wave emis- 
sion from nonconstant specific angular momentum tori undergoing 
axisymmetric oscillations while orbiting around Kerr black holes. 
The angular momentum distribution has been chosen to be increas- 
ing outward with the radial distance following a power-law de- 
pendence. The self-gravity of the discs has been neglected and 
the accretion of mass and angular momentum has been assumed 
not to affect the underlying background metric. We have also re- 
moved the possible role played by the runaway instability by se- 
lecting indeces in the power-law distribution of the specific angu- 
lar mo mentum that are above the critical value for the instabil- 
ity (cf. Daiane & Font 12004)). The results presented in this pa- 
per ha ve thus extended the previous investigation of IZanotti et alJ 
(2003) where constant specific angular momentum discs around 
Schwarzschild black holes were considered. 

A comprehensive sample of initial (marginally stable) equi- 
librium tori has been built. These tori have been subsequently per- 
turbed in various ways in order to study their dynamical response 
to the perturbations during long-term time evolutions (extending 
up to 100 orbital periods for each model). The time evolutions are 
characterized by a distinctive quasi-periodic pattern present in all 
fluid quantities as well as in the mass and angular momentum ac- 
cretion rates and in the gravitational waveforms (which have been 
extracted using the Newtonian quadrupole formula). Our study has 
been carried out using two complementary tools: an axisymmet- 
ric nonlinear hydrodynamics code and a linear perturbative code 
for vertically integrated axisymmetric tori. While the first one has 
proven useful to investigate nonlinear effects in the axisymmetric 
oscillation, the second approach has been useful in confirming and 
interpreting the numerical simulations. 

The main results of our study are as follows: The power-law 
angular momentum distribution results in more stable tori than 
in the uniform angular momentum case, yielding systematically 
smaller values of the mass accretion rate as the power-law in- 
dex increases. This is in agree ment with the p revious studies of 
iRezzolla et ail <2003bh and of iDaigne & Fond j20o4> . The com- 
parison of the eigenfrequencies calculated through the perturbative 
analysis and those obtained from the power spectra of the hydrody- 
namical variables in the numerical simulations have shown a good 
agreement between the two approaches for the first two modes / 
and 01 . They have been found close to a simple sequence of inte- 
gers 2:3, but with deviations that become significant as the power- 
law index q increases. Namely, for q = 0.2 the ratio f/01 dif- 
fers from 2/3 by ~ 15%. On the other hand, the second overtone 
02 predicted by the perturbative analysis could be detected only in 
two cases among the several power spectra calculated from the hy- 
drodynamical models. This is either due to the limitations of the 
vertically integrated assumption in the perturbative analysis in cap- 
turing higher frequency oscillations (something which is also con- 
firmed when performing selective excitations of modes), or by the 
very small power confined in these modes which prevents them to 
become visible in the power spectra. 

The comparison has also shown that the peak found at 2/ in 
the power spectra should not be interpreted as a proper eigenmode 
of the system. Rather, it represents the result of a nonlinear coupling 
of the lower-order modes as it is common in nonlinear systems sub- 
ject to nonharmonic oscillations. 

Finally, we have reported on the detectability of the oscil- 
lations of perturbed tori via their gravitational wave emission. In 
good agreement with what we already found for constant angular 
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Table 3. Computed estimates regarding the detection of the gravitational wave signals emitted by relativistic tori around Kerr black holes with power-law 
distributions of the angular momentum. From left to right the table reports the characteristic frequency, the characteristic amplitude, and the signal-to-noise 
ratio computed for three detectors, LIGO I, VIRGO, and Advanced LIGO, assuming a galactic distance for the first two, and an extragalactic distance for the 
Advanced LIGO detector. Tuf e = 100 orbital periods. 



Model 


fc (Hz) 


fc (Hz) 


fc (Hz) 


h c 


h c 


h c 


S/N 


S/N 


S/N 




LIGO I 


VIRGO 


ADV. LIGO 


LIGO I 


VIRGO 


ADV. LIGO 


LIGO I 


VIRGO 


ADV. LIGO 




(10 Kpc) 


(10 Kpc) 


(20 Mpc) 


(lOKpc) 


(lOKpc) 


(20 Mpc) 


(10 Kpc) 


(10 Kpc) 


(20 Mpc) 


Ala 


223 


236 


223 


1.5 x 10~ 20 


1.5 x 10" 20 


7.1 x 10~ 24 


29.1 


23.5 


0.26 


A2a 


267 


274 


264 


3.8 x 10- 20 


3.8 x 10" 20 


1.9 x 10~ 23 


62.6 


55.2 


0.59 


A3a 


296 


308 


291 


1.3 x 10" 20 


1.4 x 10" 20 


6.3 x 10" 24 


18.9 


17.7 


0.17 


A4a 


121 


117 


122 


1.4 x 10" 21 


1.7 x 10~ 21 


8.3 x 10" 25 


4.1 


3.2 


0.03 


A5a 


205 


207 


206 


4.6 x 10" 21 


4.7 x 10" 21 


2.3 x 10" 24 


10.1 


7.5 


0.09 


Cla 


269 


388 


237 


3.3 x 10" 22 


4.2 x 10~ 22 


1.1 x 10~ 25 


0.5 


0.5 


< 10" 2 


Dla 


446 


454 


436 


1.2 x 10" 19 


1.2 x 10" 19 


5.6 x 10" 23 


100 


122 


0.6 


D2a 


366 


393 


349 


4.3 x 10" 20 


4.5 x 10" 20 


1.8 x 10" 23 


48 


52 


0.3 


D3a 


344 


367 


336 


1.6 x 10" 20 


1.7 x 10~ 20 


7.3 x 10~ 24 


19.5 


20.4 


0.14 


D4a 


498 


505 


487 


2.3 x 10~ 19 


2.4 x 10- 19 


1.1 x 10~ 22 


175 


225 


1.0 


Ela 


272 


341 


258 


1.6 x 10" 20 


1.9 x 10" 20 


6.4 x 10" 24 


26.4 


23.8 


0.21 


E2a 


249 


323 


236 


8.5 x 10" 21 


1.1 x 10" 20 


3.3 x 10" 24 


15.1 


13.2 


0.12 


E3a 


538 


556 


506 


6.5 x 10" 20 


6.7 x 10" 20 


3.0 x 10" 23 


44 


58 


0.23 


E4a 


454 


193 


406 


2.3 x 10" 20 


2.4 x 10~ 20 


1.0 x 10" 23 


20 


24 


0.12 


E5a 


412 


191 


347 


1.1 x 10" 20 


1.3 x 10" 20 


4.3 x 10" 24 


10.9 


12.3 


0.07 



momentum discs, the chances for detection of gravitational radia- 
tion from nonconstant angular momentum tori, when sufficiently 
compact and dense, are good and within the sensitivity curves of 
LIGO and VIRGO for galactic sources, but only marginal even for 
Advanced LIGO for extragalactic sources located at 20 Mpc. 

Extensions of the work reported here include the incorporation 
of additional physics in the models such as viscosity and magnetic 
fields and will be reported in forthcoming papers. 
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APPENDIX A: DESCRIPTION OF THE ATMOSPHERE 
USED IN THE NUMERICAL SIMULATIONS 

We model the low density atmosphere surrounding the tori in terms 
of a semi-analytic solution of the relativistic stationary accretion of 
a rotating fluid onto a Kerr black hole. In this sense, it represents 
the extension of the relativistic spherical accretion solution onto a 
Schwarzschild black hole l Mich el7972t) to account for the rotation 
of both the infalling fluid and of the background spacetime. 

In practice, in addition to the continuity equation, V a (pu a ) = 
0, and to the energy conservation equation, V a T" = 0, one also 
has to ensure the conservation of the angular momentum, V a T£ = 
0. During its motion, a fluid with a general four velocity u a = 
(it*, it r , 0, tt 1 ^) follows a spiral ending at the black hole horizon 



along cones of constant 9. This motion can be derived from the 
following three constraints coming directly from the conservation 
equations 



-ghpuut = 
-ghpuu^ = 



a, 

C*2, 



(Al) 
(A2) 
(A3) 



where Ci, C2 and C3 can only depend on the polar angle 8. Di- 
viding Eqs. <A3> and <A2t gives the condition of constant specific 
angular momentum £ = —u^/ut, while the division of Eqs. iA2\ 
and JA1> provides the relativistic Bernoulli equation 

hut = C , (A4) 

where C = C2/C1. Following M ichel i 19721) . Eq. <A4t is squared 
and then differentiated to obtain, with the help of the continuity 
equation, Eq. <AH . the so-called "wind equation" 



2 du 
u 
dr 



-V 2 G(r,u) + 



H{r,u) 
[-W 2 G(r,u) + S{r,u)] = 
where it is a shorthand notation for u 



+ 
0, 

while 



(A5) 



F(r,u) 
H(r,u) 
G(r,u) 
S(r,u) 



,2 d>4> 



(A6) 
(A7) 
(A8) 
(A9) 
(A10) 



1 + g rr ii 

g u -2£g^+rg 
F{r,u)/H[r,u) 
r[u 2 d r g rr H(r, it) — F(r, it) 
{d r g u - 2£d r g t4 ' + fdrg^/H^u) 

and V 2 = dln(/ip)/dln p— 1 is the square of the sound speed. Im- 
posing the vanishing of the terms in the square brackets of Eq. <A5> 
allows to determine the position of the critical point of the flow, 
thus guaranteeing that the derivative du/dr is always finite. Note 
that according to Eq. <A5> the sound speed at the critical point is 
given by 



V 



(All) 
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which is different from the velocity of the fluid v 2 = <?;•, ?/ u J . This 
means that the critical point is not a transonic point as in the case 
of purely spherical accretion. 

In the numerical implementation of this solution we further 
impose the simplifying assumption u$ /u t = 0, thus reducing con- 
siderably the algebra. Once the density p c at the critical point has 
been given as a free parameter, the rest of the solution is computed 
as follows: 

(i) Compute V c from the equation of state. 

(ii) Compute the relation between u c and r c by equating the 
terms in the square brackets of Eq. <A5> . 

(iii) Compute the position of the critical point r c by solving 
Eq. lATTV 

(iv) Compute the constants C\,Ci, and C3 from the known val- 
ues of the solution at the critical point. 

(v) Solve the relativistic Bernoulli equation as an algebraic 
equation in the unknown it(r). Complete the solution by comput- 
ing p(r) from the continuity equation and all of the other thermo- 
dynamic quantities from the equation of state. 
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